##Tugba Bozcaga and Kelly Pasolli
##Gov 2001 Project
library(plm)
library(foreign)
library(lmtest)
library(sandwich)
library(zoo)
library(mfx)
library(xtable)
my.d<-read.dta("my.d.final.dta")

my.d$dems<-ifelse(my.d$firstpid3==1, 1, 0)
my.d$reps<-ifelse(my.d$firstpid3==2, 1, 0)

############################################################################################
###########  Probit Models with marginal effects and clustered standard errors #############
#######################Table I - HETEROGENOUS TREATMENT BY INCOME ##########################
#### Column 1 #####
probitmfx1i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob * income2 + 
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 
probitmfx1i 

##### Column 2 #####
probitmfx2i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps +
                         HH_muchlessincome * income2 + 
                         prevwelfa_ + stillunemployed + newjob + notinlabormarket +
                         logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx2i

##### Column 3 #####
probitmfx3i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lessecure * income2 + 
                         prevwelfa_ + stillunemployed + newjob + notinlabormarket + 
                         + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx3i 

##### Column 4 #####
probitmfx4i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx4i

##Column 5
probitmfx5i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + splostjob_ + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())
probitmfx5i

##Column 6
probitmfx6i<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + splostjob_ + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + unemployment + stateXwave  , data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())
probitmfx6i
############################################################################################
#######################Table II - HETEROGENOUS TREATMENT BY EDUCATION#######################
#### Column 1 #####
probitmfx1e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps +
                         lostjob * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 
probitmfx1e

##### Column 2 #####
probitmfx2e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx2e

##### Column 3 #####
probitmfx3e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lessecure * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx3e 

##### Column 4 #####
probitmfx4e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx4e

##### Column 5 #####
probitmfx5e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx5e

####Column 6 ######
probitmfx6e<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2 + unemployment + stateXwave, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx6e

############################################################################################
#######################Table III - HETEROGENOUS TREATMENT BY GENDER#########################
#### Column 1 #####

my.d$female_dummy<-ifelse(my.d$gender_=="Female", 1, 0)

probitmfx1f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob*female_dummy+
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, 
                       atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 
probitmfx1f 


##### Column 2 #####
probitmfx2f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * female_dummy +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age +  marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx2f

##### Column 3 #####
probitmfx3f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lessecure * female_dummy +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx3f 

##### Column 4 #####
probitmfx4f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : female_dummy +
                         HH_muchlessincome : female_dummy + 
                         lessecure : female_dummy +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + female_dummy + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx4f

####Column5####
probitmfx5f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : female_dummy +
                         HH_muchlessincome : female_dummy + 
                         lessecure : female_dummy + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + female_dummy + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx5f




####Column6####
probitmfx6f<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : gender_ +
                         HH_muchlessincome : gender_ + 
                         lessecure : gender_ + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2 + unemployment + stateXwave, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx6f
############################################################################################
#######################Table IV - HETEROGENOUS TREATMENT BY AGE#############################
#### Column 1 #####
probitmfx1y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob * age +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, 
                       robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 
probitmfx1y

##### Column 2 #####
probitmfx2y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * age +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx2y 

##### Column 3 #####
probitmfx3y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lessecure * age +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx3y 

##### Column 4 #####
probitmfx4y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx4y 

####Column5####
probitmfx5y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+ splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx5y

####Column6####
probitmfx6y<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+ splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2 +unemployment + stateXwave, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx6y 

############################################################################################
#######################Table V - HETEROGENOUS TREATMENT BY RACE############################
#### Column 1 #####
my.d$racedummy<-ifelse(my.d$race2=="Black"|my.d$race2=="Hispanic"|my.d$race2=="Middle Eastern"|
                         my.d$race2=="Native American"|my.d$race2=="Other", 1, 0)

probitmfx1r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lostjob * racedummy +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + age + marriage2 + 
                         HHincnotreported + wave2, data=my.d, atmean = TRUE, 
                       robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx1r

##### Column 2 #####
probitmfx2r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * racedummy + newjob +
                         prevwelfa_ + stillunemployed + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx2r

##### Column 3 #####
probitmfx3r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + 
                         lessecure * racedummy + newjob +
                         prevwelfa_ + stillunemployed + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2+ HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
                       control  = list()) 

probitmfx3r 

##### Column 4 #####
probitmfx4r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + racedummy +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx4r 

##### Column 5 #####
probitmfx5r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + racedummy+
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx5r 

##### Column 5 #####
probitmfx6r<-probitmfx(bi_welfare ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + racedummy+
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2 + unemployment +stateXwave, data=my.d, atmean = TRUE, robust = TRUE,
                       clustervar1 =   "caseid", control = list())

probitmfx6r 

############################################################################################
##################ANNEX I: Check OLS Models with  clustered standard errors ################
#######################Table I - HETEROGENOUS TREATMENT BY INCOME ##########################
#### Column 1 #####
lm1i<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob * income2 + 
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm1i)
str(lm1i)
vcov(lm1i)
lm1ic<-coeftest(lm1i,vcov=vcovHC(lm1i, cluster="my.d$caseid"), df=4309) 

##### Column 2 #####
lm2i<-lm(welfare_ ~ 
                         dems + 
                         reps +
                         HH_muchlessincome * income2 + 
                         prevwelfa_ + stillunemployed + newjob + notinlabormarket +
                         logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm2i)
lm2ic<-coeftest(lm2i,vcov=vcovHC(lm2i, cluster="my.d$caseid"), df=4342) 

##### Column 3 #####
lm3i<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lessecure * income2 + 
                         prevwelfa_ + stillunemployed + newjob + notinlabormarket + 
                         + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm3i)
lm3ic<-coeftest(lm3i,vcov=vcovHC(lm3i, cluster="my.d$caseid"), df=4342)

##### Column 4 #####
lm4i<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm4i)
lm4ic<-coeftest(lm4i,vcov=vcovHC(lm4i, cluster="my.d$caseid"), df=4305)


####Column 5####
lm5i<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + splostjob_ + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)
summary(lm5i)
lm5ic<-coeftest(lm5i,vcov=vcovHC(lm5i, cluster="my.d$caseid"), df=4304)

####Column 6####
lm6i<-lm(welfare_~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : income2 +
                         HH_muchlessincome +
                         HH_muchlessincome : income2 + 
                         lessecure +
                         lessecure : income2 + splostjob_ + income2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + unemployment + stateXwave  , data=my.d)
summary(lm6i)
lm6ic<-coeftest(lm6i,vcov=vcovHC(lm6i, cluster="my.d$caseid"), df=4232)

###############################################################################################
#######################Table II - HETEROGENOUS TREATMENT BY EDUCATION #########################
#### Column 1 #####
lm1e<-lm(welfare_ ~ 
                         dems + 
                         reps +
                         lostjob * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm1e)
lm1ec<-coeftest(lm1e,vcov=vcovHC(lm1e, cluster="my.d$caseid"), df=4557)

##### Column 2 #####
lm2e<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm2e)
lm2ec<-coeftest(lm2e,vcov=vcovHC(lm1e, cluster="my.d$caseid"), df=4592)

##### Column 3 #####
lm3e<-lm(welfare_~ 
                         dems + 
                         reps + 
                         lessecure * ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm3e)
lm3ec<-coeftest(lm3e,vcov=vcovHC(lm3e, cluster="my.d$caseid"), df=4592)

##### Column 4 #####
lm4e<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm4e)
lm4ec<-coeftest(lm4e,vcov=vcovHC(lm4e, cluster="my.d$caseid"), df=4553)

##### Column 5 #####
lm5e<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm5e)
lm5ec<-coeftest(lm5e,vcov=vcovHC(lm5e, cluster="my.d$caseid"), df=4552)

####Column 6######
lm6e<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob +
                         lostjob : ed2 +
                         HH_muchlessincome +
                         HH_muchlessincome : ed2 + 
                         lessecure +
                         lessecure : ed2 + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + gender_ + marriage2 + ed2+
                         race2 + HHincnotreported + wave2 + unemployment + stateXwave, data=my.d)

summary(lm6e)
lm6ec<-coeftest(lm6e,vcov=vcovHC(lm6e, cluster="my.d$caseid"), df=4474)

###############################################################################################
#######################Table III - HETEROGENOUS TREATMENT BY GENDER ###########################
##### Column 1 #####
lm1f<-lm(welfare_~ 
                         dems + 
                         reps + 
                         lostjob * gender_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + age + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm1f)
lm1fc<-coeftest(lm1f,vcov=vcovHC(lm1f, cluster="my.d$caseid"), df=4558)

##### Column 2 #####
lm2f<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * gender_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age +  marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm2f)
lm2fc<-coeftest(lm2f,vcov=vcovHC(lm2f, cluster="my.d$caseid"), df=4592)

##### Column 3 #####
lm3f<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lessecure * gender_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm3f)
lm3fc<-coeftest(lm3f,vcov=vcovHC(lm3f, cluster="my.d$caseid"), df=4592)

##### Column 4 #####
lm4f<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : gender_ +
                         HH_muchlessincome : gender_ + 
                         lessecure : gender_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm4f)
lm4fc<-coeftest(lm4f,vcov=vcovHC(lm4f, cluster="my.d$caseid"), df=4553)

####Column5####
lm5f<-lm(welfare_~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : gender_ +
                         HH_muchlessincome : gender_ + 
                         lessecure : gender_ + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm5f)
lm5fc<-coeftest(lm5f,vcov=vcovHC(lm5f, cluster="my.d$caseid"), df=4552)

####Column6####
lm6f<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : gender_ +
                         HH_muchlessincome : gender_ + 
                         lessecure : gender_ + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2 + unemployment + stateXwave, data=my.d)

summary(lm6f)
lm6fc<-coeftest(lm6f,vcov=vcovHC(lm6f, cluster="my.d$caseid"), df=4474)

###############################################################################################
#######################Table IV- HETEROGENOUS TREATMENT BY AGE ################################
#### Column 1 #####
lm1y<-lm(welfare_~ 
                         dems + 
                         reps + 
                         lostjob * age +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm1y)
lm1yc<-coeftest(lm6f,vcov=vcovHC(lm6f, cluster="my.d$caseid"), df=4557)

##### Column 2 #####
lm2y<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * age +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 

summary(lm2y)
lm2yc<-coeftest(lm2f,vcov=vcovHC(lm2f, cluster="my.d$caseid"), df=4592)

##### Column 3 #####
lm3y<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lessecure * age +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d) 
summary(lm3y)
lm3yc<-coeftest(lm3y,vcov=vcovHC(lm3y, cluster="my.d$caseid"), df=4592)

##### Column 4 #####
lm4y<-lm(welfare_~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)
summary(lm4y)
lm4yc<-coeftest(lm4y,vcov=vcovHC(lm4y, cluster="my.d$caseid"), df=4553)

####Column5####
lm5y<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+ splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2, data=my.d)

summary(lm5y)
lm5yc<-coeftest(lm5y,vcov=vcovHC(lm5y, cluster="my.d$caseid"), df=4552)

####Column6####
lm6y<-lm(welfare_~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : age +
                         HH_muchlessincome : age + 
                         lessecure : age+ splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2 + HHincnotreported + wave2 +unemployment + stateXwave, data=my.d)

summary(lm6y)
lm6yc<-coeftest(lm6y,vcov=vcovHC(lm6y, cluster="my.d$caseid"), df=4474)

###############################################################################################
#######################Table v - HETEROGENOUS TREATMENT BY RACE ###############################
#### Column 1 #####
my.d$racedummy<-ifelse(my.d$race2=="Black"|my.d$race2=="Hispanic"|my.d$race2=="Middle Eastern"|
                         my.d$race2=="Native American"|my.d$race2=="Other", 1, 0)

lm1r<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lostjob * racedummy +
                         prevwelfa_ + stillunemployed + newjob +
                         notinlabormarket + logusd + ed2 + age + marriage2 + 
                         HHincnotreported + wave2, data=my.d) 

summary(lm1r)
lm1rc<-coeftest(lm1r,vcov=vcovHC(lm1r, cluster="my.d$caseid"), df=4474)

##### Column 2 #####
lm2r<-lm(welfare_~ 
                         dems + 
                         reps + 
                         HH_muchlessincome * racedummy + newjob +
                         prevwelfa_ + stillunemployed + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         HHincnotreported + wave2, data=my.d) 
summary(lm2r)
lm2rc<-coeftest(lm2r,vcov=vcovHC(lm2r, cluster="my.d$caseid"), df=4564)

##### Column 3 #####
lm3r<-lm(welfare_ ~ 
                         dems + 
                         reps + 
                         lessecure * racedummy + newjob +
                         prevwelfa_ + stillunemployed + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         race2+ HHincnotreported + wave2, data=my.d) 
summary(lm3r)
lm3rc<-coeftest(lm2r,vcov=vcovHC(lm3r, cluster="my.d$caseid"), df=4592)


##### Column 4 #####
lm4r<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + racedummy +
                         prevwelfa_ + stillunemployed + newjob + 
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2, data=my.d)

summary(lm4r)
lm4rc<-coeftest(lm4r,vcov=vcovHC(lm4r, cluster="my.d$caseid"), df=4559)

##### Column 5 #####
lm5r<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + racedummy+
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2, data=my.d)

summary(lm5r)
lm5rc<-coeftest(lm5r,vcov=vcovHC(lm5r, cluster="my.d$caseid"), df=4558)

##### Column 6 #####
lm6r<-lm(welfare_ ~ 
                         dems + 
                         reps + lostjob + HH_muchlessincome + lessecure+
                         lostjob : racedummy +
                         HH_muchlessincome : racedummy + 
                         lessecure : racedummy + splostjob_ +
                         prevwelfa_ + stillunemployed + newjob + racedummy+
                         notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                         racedummy + HHincnotreported + wave2 + unemployment +stateXwave, data=my.d)
summary(lm6r)
lm6rc<-coeftest(lm6r,vcov=vcovHC(lm6r, cluster="my.d$caseid"), df=4480)

###TABLES
library(stargazer)
dont.report<-c("marriage2", "race2", "HHincnotreported", "stwa", "wave2")
stargazer(lm1ic,lm2ic,lm3ic,lm4ic,lm5ic,lm6ic, omit=dont.report,font.size="tiny")
stargazer(lm1fc,lm2fc,lm3fc,lm4fc,lm5fc,lm6fc, omit=dont.report,font.size="tiny")
stargazer(lm1ec,lm2ec,lm3ec,lm4ec,lm5ec,lm6ec, omit=dont.report,font.size="tiny")
stargazer(lm1yc,lm2yc,lm3yc,lm4yc,lm5yc,lm6yc, omit=dont.report,font.size="tiny")
stargazer(lm1rc,lm2rc,lm3rc,lm4rc,lm5rc,lm6rc, omit=dont.report,font.size="tiny")

###############################################################################################
########### Annex II - PLM Models with fixed effects and clustered standard errors ############
##### Column 1 #####
p1<-as.formula(welfare_ ~ lostjob + stillunemployed + newjob 
               + notinlabormarket + logusd + ed2 + age + logusd +
                 marriage2 + HHincnotreported + wave2)
plm1 <-plm(formula=p1, data=my.d, index=c("caseid"), model="within")
plm1c<-coeftest(plm1, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

##### Column 2 #####
p2<-as.formula(welfare_ ~ HH_muchlessincome + stillunemployed + newjob 
               + notinlabormarket +logusd + 
                 ed2 + age +  marriage2 + HHincnotreported + wave2)
plm2 <-plm(formula=p2, data=my.d, index=c("caseid"), model="within")
plm2c<-coeftest(plm2, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

##### Column 3 #####
p3<-as.formula(welfare_ ~ lessecure + stillunemployed + newjob 
               + notinlabormarket + logusd +
                 ed2 + age +  marriage2 + HHincnotreported + wave2)
plm3 <-plm(formula=p3, data=my.d, index=c("caseid"), model="within")
plm3c<-coeftest(plm3, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

##### Column 4 #####
p4<-as.formula(welfare_ ~ lostjob + HH_muchlessincome + lessecure +
                 stillunemployed + newjob + notinlabormarket + logusd + 
                 ed2 + age +  marriage2 + HHincnotreported + wave2)
plm4 <-plm(formula=p4, data=my.d, index=c("caseid"), model="within")
plm4c<-coeftest(plm4, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

##### Column 5 #####
p5<-as.formula(welfare_ ~ lostjob + HH_muchlessincome + lessecure + splostjob_ +
                 stillunemployed + newjob + notinlabormarket + logusd + 
                 ed2 + age +  marriage2 + HHincnotreported + wave2)
plm5 <-plm(formula=p5, data=my.d, index=c("caseid"), model="within")
plm5c<-coeftest(plm5, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

##### Column 6 #####
stwa<-factor(my.d$stateXwave)
p6<-as.formula(welfare_ ~ lostjob + HH_muchlessincome + lessecure + splostjob_ +
                 stillunemployed + newjob + notinlabormarket + logusd + 
                 ed2 + age +  marriage2 + HHincnotreported + wave2 + my.d$unemployment + my.d$stateXwave)
plm6 <-plm(formula=p6, data=my.d, index=c("caseid"), model="within")
summary(plm6)
plm6c<-coeftest(plm6, vcov=function(x) vcovHC(x, cluster="group", type="HC1"))

######## CREATE TABLE 1 USING STARGAZER #######
dont.report<-c("marriage2", "race2", "HHincnotreported", "stwa", "wave2")
# exclude variables that are not reported and we don't want in our table
library(stargazer)
stargazer(plm1c, plm2c, plm3c, plm4c, plm5c, plm6c, omit=dont.report)

######## RANDOM OR FIXED EFFECTS? #######
######## HAUSMAN TEST #######
phtest(plm4, plm4random)

####################### GRAPHS #######################


#### Income and Job Loss ####

glm1<- glm(bi_welfare ~ 
             dems + 
             reps + 
             lostjob * income2 + 
             prevwelfa_ + stillunemployed + newjob + 
             notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
             race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
           control  = list()) 

glm1
cov.mat<-vcov(glm1)
cov.mat
mod_frame=model.frame(glm1)
beta1<--.167
beta3<-.137
income.vals<-c(seq(2,6,0.01))
effects = beta1 + beta3*income.vals
var<-cov.mat[4,4] + (income.vals^2)*cov.mat[28,28] + 2*(income.vals)*cov.mat[4,28]
se<-sqrt(var)
z_score = qnorm(1 - ((1 - .95)/2))
upper=effects+z_score*se
lower=effects-z_score*se


plot(income.vals, effects, type="l", main="Figure 2: Relationship between Job Loss and
     Income on Support for Social Policy", xlab="Logged Household Income", ylab="Marginal Effect of Job Loss",
     cex.main=0.8, cex.lab=0.7, ylim=c(0, 0.6))
lines(y=upper, x=income.vals, lty=2)
lines(y=lower, x=income.vals, lty=2)

legend("topright",
       c("95% confidence intervals"),
       lty=c(1, 4), lwd=2, ncol=1, cex=0.5)

###### Income and Less Secure

glm2<- glm(bi_welfare ~ 
             dems + 
             reps +
             HH_muchlessincome * income2 + 
             prevwelfa_ + stillunemployed + newjob + notinlabormarket +
             logusd + ed2 + age + gender_ + marriage2 + 
             race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
           control  = list()) 

glm2
cov.mat<-vcov(glm2)
cov.mat
mod_frame=model.frame(glm2)
beta1<-.521
beta3<--.146
income.vals<-c(seq(2,6,0.01))
effects = beta1 + beta3*income.vals
var<-cov.mat[5,5] + (income.vals^2)*cov.mat[28,28] + 2*(income.vals)*cov.mat[5,28]
se<-sqrt(var)
z_score = qnorm(1 - ((1 - .95)/2))
upper=effects+z_score*se
lower=effects-z_score*se


plot(income.vals, effects, type="l", ylim=c(-0.7, 0.7), main="Figure 2: Relationship between Feeling Less
     Secure and Income on Support for Social Policy", xlab="Logged Household Income", ylab="Marginal Effect of Feeling Less Secure",
     cex.main=0.8, cex.lab=0.7)
lines(y=upper, x=income.vals, lty=2)
lines(y=lower, x=income.vals, lty=2)

legend("topright",
       c("95% confidence intervals"),
       lty=c(1, 4), lwd=2, ncol=1, cex=0.5)

